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The dwarf spheroidal satellite galaxies (dSphs) of the Milky Way are some of the most dark matter 
(DM) dominated objects known. We report on gamma-ray observations of Milky Way dSphs based 
on 6 years of Fermi Large Area Telescope data processed with the new Pass 8 event-level analysis. 

None of the dSphs are significantly detected in gamma rays, and we present upper limits on the DM 
annihilation cross section from a combined analysis of 15 dSphs. These constraints are among the 
strongest and most robust to date and lie below the canonical thermal relic cross section for DM of 
mass ($100 GeV annihilating via quark and r-lepton channels. 

PACS numbers: 95.35.+d, 95.85.Pw, 98.52.Wz 
Keywords: dark matter; gamma rays; dwarf galaxies 

INTRODUCTION are an attractive candidate to constitute some or all of 

DM [2-4]. The relic abundance of WIMPs is deter¬ 
mined by their annihilation cross section at freeze-out 
Approximately 26% of the energy density of the uni- [5], and the characteristic weak-scale cross sections of 

verse is composed of non-baryonic cold dark matter WIMPs can naturally produce a relic abundance equal 

(DM) [1], Weakly interacting massive particles (WIMPs) 
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to the observed abundance of DM. Self-annihilation of 
WIMPs would continue today in regions of high DM den¬ 
sity and result in the production of energetic Standard 
Model particles. The large mass of the WIMP (ttidm) 
permits the production of gamma rays observable by the 
Fermi Large Area Telescope (LAT), which is sensitive to 
energies ranging from 20 MeV to > 300 GeV. 

Kinematic data indicate that the dwarf spheroidal 
satellite galaxies (dSphs) of the Milky Way contain a 
substantial DM component [6, 7]. The gamma-ray signal 
flux at the LAT, (j) s (phcm _2 s _1 ), expected from the 
annihilation of DM with a density distribution pdm(^) is 
given by 


4> s { AD) 
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Here, the first term is dependent on the particle physics 
properties i.e., the thermally-averaged annihilation 
cross section, (av), the particle mass, ?tidm> and the 
differential gamma-ray yield per annihilation, cW 7 /dA 7 , 
integrated over the experimental energy range. 1 The 
second term, known as the J-factor, is the line-of-sight 
(l.o.s.) integral through the DM distribution integrated 
over a solid angle, Afi. 

Milky Way dSphs can give rise to J-factors in excess of 
10 19 GeV 2 cm -5 [8, 9], which, coupled with their lack of 
non-thermal astrophysical processes, makes them good 
targets for DM searches via gamma rays. Gamma-ray 
searches for dSphs yield some of the most stringent con¬ 
straints on (av), particularly when multiple dSphs are 
analyzed together using a joint likelihood technique [10- 
15]. Limits on (av) derived from observations of dSphs 
have begun to probe the Iow-TOdm parameter space for 
which the WIMP abundance matches the observed DM 
relic density. 

In contrast, DM searches in the Galactic center take 
advantage of a J-factor that is 0(100) times larger, al¬ 
though gamma-ray emission from non-thermal processes 
makes a bright, structured background. Several stud¬ 
ies of the Galactic center interpret an excess of gamma 
rays with respect to modeled astrophysical backgrounds 
as a signal of 20 to 50 GeV WIMPs annihilating via the 
bb channel [16-19]. Coincidentally, the largest deviation 


1 Strictly speaking, the differential yield per annihilation in Equa¬ 
tion (1) is a sum of differential yields into specific final states: 
dTV 7 /d-E 7 = Bj d nJ) /di? 7 , where Bf is the branching frac¬ 
tion into final state /. Here, we make use of Equation (1) in the 
context of single final states only. 


from expected background in some previous studies of 
dSphs occurred for a similar set of WIMP characteris¬ 
tics; however, this deviation was not statistically signifi¬ 
cant [13]. 

Using a new LAT event-level analysis, known as 
Pass 8, we re-examine the sample of 25 Milky Way 
dSphs from Ackermann et al. [13] using six years of 
LAT data. The Pass 8 data benefits from an improved 
point-spread function (PSF), effective area, and energy 
reach. More accurate Monte Carlo simulations of the 
detector and the environment in low-Earth orbit have 
reduced the systematic uncertainty in the LAT instru¬ 
ment response functions (IRFs) [20]. Within the stan¬ 
dard photon classes, Pass 8 offers event types, subdivi¬ 
sions based on event-by-event uncertainties in the direc¬ 
tional and energy measurements, which can increase the 
sensitivity of likelihood-based analyses. In this work we 
use a set of four PSF event-type selections that sub¬ 
divide the events in our data sample according to the 
quality of their directional reconstruction. In addition 
to the improvements from Pass 8, we employ the up¬ 
dated third LAT source catalog (3FGL), based on four 
years of Pass 7 Reprocessed data, to model point-like 
background sources [21]. Together, these improvements, 
along with an additional two years of data taking, lead 
to a predicted increase in sensitivity of 70% relative to 
the four-year analysis of Ackermann et al. [13] for the bb 
channel at 100 GeV. More details on Pass 8 and other 
aspects of this analysis can be found in Supplemental 
Material [22]. 


LAT DATA SELECTION 

We examine six years of LAT data (2008-08-04 to 2014- 
08-05) selecting Pass 8 SOURCE-class events in the en¬ 
ergy range between 500 MeV and 500 GeV. We selected 
the 500 MeV lower limit to mitigate the impact of leakage 
from the bright limb of the Earth because the PSF broad¬ 
ens considerably below that energy. To further avoid 
contamination from terrestrial gamma rays, events with 
zenith angles larger than 100° are rejected. We also re¬ 
move time intervals around bright GRBs and solar flares 
following the prescription used for the 3FGL catalog. We 
extract from this data set 10° x 10° square regions of in¬ 
terest (ROIs) in Galactic coordinates centered at the po¬ 
sition of each dSph specified in Table I. 

At a given energy, 20%-40% of the events classified as 
photons in our six-year Pass 8 data set are shared with 
the analysis of Ackermann et al. [13]. The low fraction of 
shared events can be attributed primarily to the larger 
time range used for the present analysis (four versus six 
years) and the increase in gamma-ray acceptance of the 
P8R2 .SOURCE event class relative to P7REP .CLEAN. 
At most, the Pass 7 events can represent 35%-50% of the 
new, larger sample. Migration of the individual recon- 
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structed events, particularly residual cosmic rays, across 
ROI and class selection boundaries further reduces the 
overlap, making the two analyses nearly statistically in¬ 
dependent [22]. 


J-FACTORS FOR DWARF SPHEROIDAL 
GALAXIES 


The DM content of dSphs can be determined through 
dynamical modeling of their stellar density and velocity 
dispersion profiles [23-25]. Recent studies have shown 
that an accurate estimate of the dynamical mass of 
a dSph can be derived from measurements of the av¬ 
erage stellar velocity dispersion and half-light radius 
alone [26, 27]. The total mass within the half-light radius 
and the integrated J-factor have been found to be fairly 
insensitive to the assumed DM density profile [13, 25, 28]. 
We assume that the DM distribution in dSphs follows a 
Navarro-Frenk-White (NFW) profile [29], 


Pdm(V) 


Por 3 s 

r(r s + r) 


2 ’ 


( 2 ) 


where r s and po are the NFW scale radius and charac¬ 
teristic density, respectively. We take J-factors and other 
physical properties for the Milky Way dSphs from Ack- 
ermann et al. [13] (and references therein). 


DATA ANALYSIS 

We perform a binned Poisson maximum-likelihood 
analysis in 24 bins of energy, 2 logarithmically spaced 
from 500 MeV to 500 GeV, and an 0.1° angular pixeliza- 
tion. The low-energy bound of 500 MeV is selected to 
mitigate the impact of leakage from the bright limb of 
the Earth because the PSF broadens considerably be¬ 
low that energy. The high-energy bound of 500 GeV is 
chosen to mitigate the effect of the increasing residual 
charged-particle background at higher energies [30]. The 
data were analyzed with the Fermi Science Tools 3 ver¬ 
sion 10-01-01 and the P8R2_SOURCE_V6 IRFs. Our 
diffuse background model includes a structured Galac¬ 
tic component and a spatially isotropic component that 
represents both extragalactic emission and residual par¬ 
ticle contamination. 4 Because the energy resolution of 
the LAT was not accounted for when fitting the Galac¬ 
tic diffuse model, differences in response (energy resolu¬ 
tion and effective area) between IRF sets lead to different 


2 Constraints are insensitive to finer binning. 

2 http://fermi.gsfc.nasa.gov/ssc/data/analysis/software 

4 http://fermi.gsfc.nasa.gov/ssc/data/access/lat/ 
BackgroundModels.html 


measured intensities for this component. Thus, a small 
energy-dependent scaling has been applied to the Pass 7 
Reprocessed Galactic diffuse model. Changes with re¬ 
spect to the Pass 7 Reprocessed model are less than 5% 
above 100 MeV. Details on the derivation of the rescaled 
model are given in [22]. The gamma-ray characteristics 
of nearby point-like sources are taken from the 3FGL 
catalog [21], 

We perform a bin-by-bin likelihood analysis of the 
gamma-ray emission coincident with each dSph follow¬ 
ing the procedure of Ackermann et al. [13]. The flux 
normalizations of the Galactic diffuse and isotropic com¬ 
ponents and 3FGL catalog sources within the 10° x 10° 
ROI were fit simultaneously in a binned likelihood anal¬ 
ysis over the broadband energy range from 500 MeV to 
500 GeV. The normalizations of the background sources 
are insensitive to the inclusion of a putative power-law 
source at the locations of the dSphs, which is consis¬ 
tent with the lack of any strong signal associated with 
the dSphs. Fixing the normalizations of the background 
sources with the broad-band fit before fitting each bin in¬ 
dividually avoids numerical instability resulting from the 
fine binning in energy and the degeneracy of the diffuse 
background components at high Galactic latitudes. 

After fixing the background normalizations, we scan 
the likelihood as a function of the flux normalization of 
the putative DM signal independently in each energy bin 
(this procedure is similar to that used to evaluate the 
spectral energy distribution of a source). Within each 
bin, we model the putative dSph source with a power- 
law spectral model (dN/dE oc E~ r ) with spectral index 
of r = 2. By analyzing each energy bin separately, we 
avoid selecting a single spectral shape to span the entire 
energy range at the expense of introducing additional 
degrees of freedom into the fit. 

While the bin-by-bin likelihood function is essentially 
independent of spectral assumptions, it does depend on 
the spatial model of the DM distribution in the dSphs. 
We model the dSphs with spatially extended NFW DM 
density profiles projected along the line of sight. The 
angular extent of the emission profile for each dSph is 
set by the scale radius of its DM halo, which contains 
approximately 90% of the total annihilation flux. We 
use the set of DM halo scale radii from Ackermann et al. 
[13], which span a range of subtended angles between 0.1° 
and 0.4°. 

We test a wide range of DM annihilation hypotheses 
by using predicted gamma-ray spectra to tie the signal 
normalization across the energy bins. Spectra for DM 
annihilation are generated with the DMFIT package based 
on Pythia 8.165 [13, 31, 32], We reconstruct a broad¬ 
band likelihood function by multiplying the bin-by-bin 
likelihood functions evaluated at the predicted fluxes for 
a given DM model. 

We combine the broad-band likelihood functions across 
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15 of the observed dSphs 5 and include statistical uncer¬ 
tainties on the J-factors of each dSph by adding an ad¬ 
ditional J-factor likelihood term to the binned Poisson 
likelihood for the LAT data. The J-factor likelihood for 
target i is given by 

Wi I ^obs,ijOi) = “ . . 7= 

m(10) Jobs.iV 27T(Tj (3) 

^ g—(l°glO W>) — loglO C^ob8,i)) 2 /2CT^ 

where Ji is the true value of the J-factor and J 0 bs,i is 
the measured J-factor with error or This parameteriza¬ 
tion of the J-factor likelihood is obtained by fitting a log¬ 
normal function with peak value J 0 bs,i to the posterior 
distribution for each J-factor as derived by Martinez [8], 
providing a reasonable way to quantify the uncertainties 
on the J-factors. This approach is a slight modification 
of the approach in Ackermann et al. [10, 13], where an 
effective likelihood was derived considering a flat prior 
on the J-factors. We note that the J-factor correction is 
only intended to incorporate the statistical uncertainty 
in the J-factors, and not the systematic uncertainty re¬ 
sulting from the fitting procedure or choice of priors [22]. 
More details on the derivation of the J-factor likelihood 
and the effects of systematic uncertainties can be found 
in Supplemental Material [22]. 

Combining the broad-band gamma-ray and J-factor 
likelihood functions, our likelihood function for target i 
becomes, 

A(/Li,0i = {aj, Ji} \T>i) =C i {n,e i I T>i)Cj{Ji I J 0 bs,i> &i). 

(4) 

Here, fi are the parameters of the DM model, 0i is the 
set of nuisance parameters that includes both nuisance 
parameters from the LAT analysis (a:;) and the dSph J- 
factor (A), and Vi is the gamma-ray data. We incorpo¬ 
rate additional information about the event-wise quality 
of the angular reconstruction by forming the LAT likeli¬ 
hood function (£j) from the product of likelihood func¬ 
tions for four PSF event types. The four PSF event types 
(PSFO, PSF1, PSF2, and PSF3) subdivide the events 
in the SOURCE-class data set into exclusive partitions 
(Vi j) in order of decreasing uncertainty on the direction 
measurement. The resulting joint LAT likelihood func¬ 
tion is given by 

A(M) Oi I Vi) = Ci(i Li, 0, I V itj ). 

j 

The spectral and spatial model of gamma-ray counts for 
each event type partition is evaluated using a set of IRFs 
computed for that class and type selection. 


We evaluate the significance of DM hypotheses using 
a test statistic (TS) defined as 


TS = -2 In 


( C{^.,6\V) \ 

\c{jx,e\v)) ’ 


( 6 ) 


where /x 0 are the parameters of the null (no DM) hy¬ 
pothesis and jj, and 6 are the best-fit parameters un¬ 
der the DM hypothesis. £ can here be either the like¬ 
lihood for an individual dSph or the joint likelihood for 
the dSphs in our combined sample. We note that fol¬ 
lowing the methodology of Ackermann et al. [13] we use 
background parameters (6) derived under the hypothe¬ 
sis of a DM source with a T = 2 power-law spectrum 
when evaluating both the null and DM hypotheses. This 
is a good approximation as long as the best-fit signal is 
small relative to the background in the ROI. Based on 
the asymptotic theorem of Chernoff [33], the TS can be 
converted to a significance based on a mixture of y 2 dis¬ 
tributions. The validity of this assumption is examined 
further in Supplemental Material [22]. 


RESULTS 

We find no significant gamma-ray excess associated 
with the Milky Way dSphs when analyzed individually 
or as a population. In the combined analysis of 15 dSphs, 
the largest deviation from the background-only hypothe¬ 
sis has TS = 1.3 occurring for todm = 2 GeV annihilating 
through the e + e _ channel. Among the dSphs in our com¬ 
bined analysis, the dSph with the largest individual sig¬ 
nificance is Sculptor with TS = 4.3 for todm = 5 GeV an¬ 
nihilating through the fi~ channel. The maximum TS 
of our combined analysis is well below the threshold set 
for gamma-ray source detection and is completely con¬ 
sistent with a background fluctuation [21]. We set upper 
limits on (av) at 95% confidence level (CL) for WIMPs 
with todm between 2 GeV and 10 TeV annihilating into 
six different standard model channels (66, r + r _ , 
e + e _ , W + W ~, uu). 6 Figure 1 shows the comparison of 
the limits for the 66 and t + t~ channels with expectation 
bands derived from the analysis of 300 randomly selected 
sets of blank fields. Sets of blank fields are generated by 
choosing random sky positions with |6| > 30° that are 
centered at least 0.5° from 3FGL catalog sources. We ad¬ 
ditionally require fields within each set to be separated by 
at least 7°. Our expected limit bands are evaluated with 
the 3FGL source catalog based on four years of Pass 7 
Reprocessed data and account for the influence of new 
sources present in the six-year Pass 8 data set. 


5 Selected to have kinematically determined J-factors and avoid 
ROI overlap. The set is identical to that in Ackermann et al. 

[13]. 


6 Results for all channels as well as bin-by-bin likelihood func¬ 
tions for each target are available in machine-readable format 
at: http://www-glast.stanford.edu/pub_data/1048/. 
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Comparing with the results of Ackermann et al. [13], 
we find a factor of 3-5 improvement in the limits for all 
channels using six years of Pass 8 data and the same 
sample of 15 dSphs. The larger data set as well as the 
gains in the LAT instrument performance enabled by 
Pass 8 both contribute to the increased sensitivity of the 
present analysis. An additional 30-40% improvement in 
the limit can be attributed to the modified functional 
form chosen for the J-factor likelihood (Equation 3). Sta¬ 
tistical fluctuations in the Pass 8 data set also play a 
substantial role. Because the Pass 8 six-year and Pass 7 
Reprocessed four-year event samples have a shared frac¬ 
tion of only 20-40%, the two analyses are nearly statisti¬ 
cally independent. For masses below 100 GeV, the upper 
limits of Ackermann et al. [13] were near the 95% upper 
bound of the expected sensitivity band while the limits 
in the present analysis are within one standard deviation 
of the median expectation value. 

Uncertainties in the LAT IRFs, modeling of the diffuse 
background, and estimation of J-factors all contribute 
systematic errors to this analysis. By examining maximal 
variations of each contributor, we find that at 100 GeV 
they lead to ±9%, ±8%, and ±33% shifts in our limits, 
respectively (see Supplemental Material [22]). 

Our results begin to constrain some of the preferred pa¬ 
rameter space for a DM interpretation of a gamma-ray 
excess in the Galactic center region [16-19]. As shown 
in Figure 2, for interpretations assuming a bb final state, 
the best-fit models lie in a region of parameter space 
slightly above the 95% CL upper limit from this anal¬ 
ysis, with an annihilation cross section in the range of 
(l-3)xl0~ 26 cm 3 s -1 and TOdm between 25 and 50 GeV. 
However, uncertainties in the structure of the Galactic 
DM distribution can significantly enlarge the best-fit re¬ 
gions of (av), channel, and todm [38]. 

In conclusion, we present a combined analysis of 15 
Milky Way dSphs using a new and improved LAT data 
set processed with the Pass 8 event-level analysis. We ex¬ 
clude the thermal relic annihilation cross section (~ 2.2 x 
10~ 26 cm 3 s _1 ) for WIMPs with todm 100 GeV annihi¬ 
lating through the quark and r-lepton channels. Our 
results also constrain DM particles with todm above 
100 GeV surpassing the best limits from Imaging Atmo¬ 
spheric Cherenkov Telescopes for masses up to ~ 1 TeV 
for quark channels and ~ 300 GeV for the r-lepton chan¬ 
nel. These constraints include the statistical uncertainty 
on the DM content of the dSphs. The future sensitivity 
to DM annihilation in dSphs will benefit from additional 
LAT data taking and the discovery of new dSphs with 
upcoming optical surveys such as the Dark Energy Sur¬ 
vey [39] and the Large Synoptic Survey Telescope [40]. 


TABLE I. Properties of Milky Way dSphs. 


Name 


fe a 

Distance 

Jogio 

(j 

obs) b 

Ref. 


(deg) 

(deg) 

(kpc) 

(l°Slo[ GeV 2 cm 

■ 5 1) 

Bootes I 

358.1 

69.6 

66 

18.8 

± 

0.22 

[41] 

Canes Venatici II 

113.6 

82.7 

160 

17.9 

± 

0.25 

[42] 

Carina 

260.1 

-22.2 

105 

18.1 

± 

0.23 

[43] 

Coma Berenices 

241.9 

83.6 

44 

19.0 

± 

0.25 

[42] 

Draco 

86.4 

34.7 

76 

18.8 

± 

0.16 

[44] 

Fornax 

237.1 

-65.7 

147 

18.2 

± 

0.21 

[43] 

Hercules 

28.7 

36.9 

132 

18.1 

± 

0.25 

[42] 

Leo II 

220.2 

67.2 

233 

17.6 

± 

0.18 

[45] 

Leo IV 

265.4 

56.5 

154 

17.9 

± 

0.28 

[42] 

Sculptor 

287.5 

-83.2 

86 

18.6 

± 

0.18 

[43] 

Segue 1 

220.5 

50.4 

23 

19.5 

± 

0.29 

[46] 

Sextans 

243.5 

42.3 

86 

18.4 

± 

0.27 

[43] 

Ursa Major II 

152.5 

37.4 

32 

19.3 

± 

0.28 

[42] 

Ursa Minor 

105.0 

44.8 

76 

18.8 

± 

0.19 

[44] 

Willman 1 

158.6 

56.8 

38 

19.1 

± 

0.31 

[47] 

Bootes II c 

353.7 

68.9 

42 


- 


- 

Bootes III 

35.4 

75.4 

47 


- 


- 

Canes Venatici I 

74.3 

79.8 

218 

17.7 

± 

0.26 

[42] 

Canis Major 

240.0 

-8.0 

7 


- 


- 

Leo I 

226.0 

49.1 

254 

17.7 

± 

0.18 

[48] 

Leo V 

261.9 

58.5 

178 


- 


- 

Pisces II 

79.2 

-47.1 

182 


- 


- 

Sagittarius 

5.6 

-14.2 

26 


- 


- 

Segue 2 

149.4 

-38.1 

35 


- 


- 

Ursa Major I 

159.4 

54.4 

97 

18.3 

± 

0.24 

[42] 


a Galactic longitude and latitude. 

b J-factors are calculated assuming an NFW density profile and 
integrated over a circular region with a solid angle of 
AQ ~ 2.4 x 10 -4 sr (angular radius of 0.5°). 
c dSphs below the horizontal line are not included in the 
combined analysis. 
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FIG. 1. Constraints on the DM annihilation cross section at 95% CL for the bb (left) and t + t~ (right) channels derived from 
a combined analysis of 15 dSphs. Bands for the expected sensitivity are calculated by repeating the same analysis on 300 
randomly selected sets of high-Galactic-latitude blank fields in the LAT data. The dashed line shows the median expected 
sensitivity while the bands represent the 68% and 95% quantiles. For each set of random locations, nominal J-factors are 
randomized in accord with their measurement uncertainties. The solid blue curve shows the limits derived from a previous 
analysis of four years of Pass 7 Reprocessed data and the same sample of 15 dSphs [13]. The dashed gray curve in this and 
subsequent figures corresponds to the thermal relic cross section from Steigman et al. [5]. 
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FIG. 2. Comparison of constraints on the DM annihilation cross section for the bb (left) and r + r“ (right) channels from this 
work with previously published constraints from LAT analysis of the Milky Way halo (3er limit) [34], 112 hours of observations 
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channel limits for the Galactic Center H.E.S.S. observations are taken from Abazajian and Harding [37] and assume an Einasto 
Milky Way density profile with pg = 0.389 GeV cm~% Closed contours and the marker with error bars show the best-fit cross 
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Supplemental Material: Searching for Dark Matter Annihilation from Milky Way 
Dwarf Spheroidal Galaxies with Six Years of Fermi-LAT Data 


Pass 8 Event-Level Analysis 

Pass 8 is a new event-level analysis for the LAT instrument and is the successor to the Pass 7 Reprocessed event- 
level analysis [20, 49, 50]. Some of the key features of Pass 8 are new algorithms to identify out-of-time signals, a new 
tree-based pattern recognition for the tracker subsystem, and an improved energy reconstruction that extends the 
LAT energy range below 100 MeV and above 1 TeV. Pass 8 implements a new classification analysis based on boosted 
decision trees (BDTs), which provides enhanced background rejection power relative to Pass 7 Reprocessed [20]. 
The Pass 8 event-level analysis enhances the capabilities of the LAT in all metrics relevant for high-level science 
analysis. In the energy range between 1 GeV and 10 GeV, the new P8R2_SOURCE event class has a 30-40% better 
point-source sensitivity than the P7REP_CLEAN event class. 

Pass 8 introduces an event type classification scheme that partitions events within a class according to their recon¬ 
struction quality. The event type classification is a generalization of the existing conversion type designation that 
identifies events converting in the Front or Back section of the tracker [51, 60]. Pass 8 defines eight new event type 
selections based on a sequence of energy-dependent cuts on the BDT variables that categorize the quality of the di¬ 
rection and energy reconstruction. Four event types categorize the quality of the directional reconstruction (PSF0 to 
PSF3), and another four do so for the energy reconstruction (EDISP0 to EDISP3). By construction, these selections 
partition the gamma-ray acceptance at each energy such that an event class will have approximately the same number 
of events of each type. 

Our maximum likelihood analysis of the dSphs combines the four P8R2_SOURCE_V6 PSF event types in a joint 
likelihood function. Although each event type contains approximately the same fraction of the total instrument 
acceptance, the angular resolution as measured by the 68% and 95% containment radii of the PSF is significantly 
better for events belonging to the best PSF event types. At 3.16 GeV the 68% (95%) containment radii of the 
acceptance-weighted PSF for the best and worst PSF event types (PSF3 and PSF0) is 0.17° (0.35°) and 0.92° (2.3°). 
By combining the event types in a joint likelihood function, we weight the contribution of events within a class by 
their reconstruction quality, e.g., events with the least well-characterized direction (PSF0) are assigned the lowest 
weight when testing the hypothesis of a putative DM source. We estimate that splitting the event sample by event 
type improves the sensitivity to an isolated point source by 10%. We expect that larger sensitivity gains will occur in 
regions where the gamma-ray intensity is strongly non-uniform. 

Given these improvements, along with two additional years of data, our flux constraints are expected to improve 
by a factor of ~1.7 below 10 GeV and ~2.2 above 100 GeV relative to the analysis of Ackermann et al. [13]. Although 
both the Pass 7 Reprocessed and Pass 8 analyses yield limits on the DM annihilation cross section within their 
respective 95% sensitivity bands, their constraints differ by a factor that is appreciably larger than expected from the 
median experimental sensitivities. For the bb channel, the Pass 8 constraints are ~ 5 times lower at 100 GeV. For two 
independent data sets, statistical fluctuations can easily account for the difference in limit realizations. 

We find that, at a given energy, only 20-40% of events in the six-year Pass 8 SOURCE-class data set are shared 
with the four-year Pass 7 Reprocessed CLEAN-class data set. If the Pass 8 SOURCE-class selection retains all 
events in the Pass 7 Reprocessed data set, we expect this fraction to equal the product of the ratio of gamma-ray 
acceptances for the two event classes with the ratio of observation times - in this case 35-50%. Since the basic event 
reconstruction in Pass 8 is fundamentally different, however, the characteristics of individual events change slightly. 
Events near the threshold of any of our analysis cuts can migrate out of our data selection and be excluded from the 
Pass 8 analysis. Specifically, 1 -3% of the Pass 7 Reprocessed events were reconstructed outside of our ROIs, and 
10-15% were not deemed likely enough to be photons to be included in the Pass 8 SOURCE class. Figure 3 shows the 
fraction of shared events as a function of energy for the set of dSphs in our combined analysis, along with a sample 
of events from the Earth limb, selected from time periods when the magnitude of the rocking angle of the LAT was 
greater than 52°. The Earth’s limb is an extremely pure photon source, and the fact that its shared fraction lies 
only slightly below the event class acceptance ratio indicates that the P8R2_SOURCE selection has a high efficiency 
for retaining gamma-ray events in the P7REP_CLEAN selection. The lower shared fraction observed in the dSph 
ROIs can be attributed to misclassified charged-particle events that constitute ~ 30-40% of the diffuse, high-latitude 
background between 1 and 10 GeV. These events typically lie near the boundary of the selections used to discriminate 
gamma rays from charged-particle backgrounds and are much more likely to migrate out of a given event class than 
a true gamma-ray event. 

Using the fractional overlap, we can estimate the evolution of a background fluctuation between analyses. For 
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FIG. 3. Fraction of events in the 6-year Pass 8 SOURCE data set that are also in the 4-year Pass 7 Reprocessed CLEAN 
data set. The lower points (black) show the shared fraction for events within 10° of one of the 15 dSphs used in the combined 
analysis. The upper points (blue) show the expectation for the shared fraction of a pure gamma-ray sample as derived from 
events taken from the Earth Limb. The red, solid line shows the acceptance ratio of P7REP_CLEAN_V15 to P8R2_SOURCE_V6 
and represents the maximum possible shared event fraction. Both the Earth Limb fraction and acceptance ratios have been 
scaled by the ratio of the observation times of the 4- and 6-year analyses. 


example, the largest Pass 7 Reprocessed excess occurred for the bb channel at masses between 10 and 25 GeV, and 
had a local significance of ~ 3<r. If we assume this excess resulted from an upward fluctuation of the background, the 
addition of new data is likely to reduce its significance. If the excess is due to a statistical fluctuation, we quantitatively 
expect the original significance, 04 , to drop by a factor that depends only on the intrinsic fraction of shared events, 
fi, and the ratio of observation times t\/t 2 , 


(1) 

where the intrinsic fraction is related to the observed fraction (/ ohs ) by fi ~ fobs(h/ti)- In the energy range between 
1 and 10 GeV, the Pass 8 analysis shares f 0 b s ~35% of its events with Pass 7 Reprocessed and has a 50% longer 
observation period corresponding to an intrinsic shared fraction of fi ~52%. A 3a background fluctuation should 
therefore drop to ~1.3c, which is consistent with our observations. 


Galactic Diffuse Background Model 

We model diffuse backgrounds in the dSph ROIs using templates that account for both astrophysical backgrounds 
and residual particle contamination. Our model for the Galactic diffuse emission is based on the Pass 7 Reprocessed 
diffuse emission model (gll_iem_v05_revl. f it) that was derived from a fit to Pass 7 Reprocessed data with the 
P7REP_CLEAN_V15 IRFs. Although models of the Galactic diffuse emission are formally independent of the IRFs, 
the Pass 7 Reprocessed model was fit without accounting for energy dispersion which introduces some dependence 
on the Pass 7 Reprocessed IRFs. To create a model that can be used self-consistently with Pass 8 data, we have 
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derived a rescaled model (gll_iem_v06. f it) that accounts for the influence of energy dispersion by correcting for the 
differences in apparent intensity in the Pass 7 Reprocessed and Pass 8 data sets. 

The correction for energy dispersion is derived from the ratio of the counts distributions computed with and without 
the correction for energy dispersion. Given an intensity, 1(E), the counts densities evaluated with energy dispersion 
(Ci) and assuming perfect energy resolution (C 2 ) are given by 

Ci(E') = J J I(E)A eS (E,9)t ohs (9)D(E’-E,9)dEd9, (2) 

C 2 (E') = 1(E) j j A cS (E, 6)t obs (e)6(E' - E)dEd9, (3) 

where E and E' are the true and measured energies, 9 is the LAT incident angle {9 = 0 is normal to the top of 
the LAT), A e g(E, 9) is the effective area, t obs is the integrated livetime, and D(E';E,9) is the energy dispersion. 
Defining ratios between the counts densities as Rp$ = C[ 8 /C P8 and Rpirep = C P7REP /C P7REP , the rescaled 
model is Ipg(E) = (Rpg/ Rp 7 rep)Ip 7 REp(E). We note that the correction depends on the true intensity (1(E)) in 
Equations 2 and 3. We find that the correction is not strongly dependent on the assumed spectrum, and we use here 
the average all-sky intensity of the Pass 7 Reprocessed model. Systematic uncertainties associated with our model 
for the Galactic diffuse emission are discussed in more detail in the Systematic Uncertainties section. 


Statistical Methodology 

We use a maximum likelihood-based statistical formalism [52] to test the DM signal hypothesis and derive confidence 
intervals on ( av ). Our global likelihood function for (av) is constructed from the product of likelihood functions for 
individual dSphs in our sample. We compute the profile likelihood function as a function of (av) by maximizing the 
global likelihood function with respect to the nuisance parameters for each dSph (6, = {aq, J;}): 

X((av),m DM ) (4) 

i 

Here, a t are the best-fit parameters derived from a global fit with a free-normalization DM source with a T = 2 
power-law spectrum and Ji is the J-factor that maximizes C t for a given (av) and todm- Confidence intervals on (av) 
for a given todm are calculated with the delta-log-likelihood technique, requiring a change in the profile log-likelihood 
of 2.71/2 from its maximum for a 95% CL upper limit [53]. 

Among the nuisance parameters in Equation (4), we distinguish between parameters constrained by the gamma-ray 
data (a) and the J-factors (J), which are constrained by an independent analysis of stellar kinematics. We use J- 
factors derived from a two-level Bayesian hierarchical modeling analysis that incorporates information on both stellar 
kinematics and priors on the distribution of global dSph properties [ 8 ]. For each dSph the Bayesian analysis provides 
a posterior distribution function, V(Ji), which we approximate with a log-normal distribution with central value, 
J 0 b s ,i, and uncertainty, a,. 

Following the approach developed in Ackermann et al. [10, 13], we account for statistical uncertainty on the J-factor 
by multiplying the LAT likelihood function with a J-factor likelihood function, Cj(Ji | </ o bs ,0 We construct an 
ansatz for the J-factor likelihood function by equating the sampling distribution of J 0 bs,i with V(Ji). With this 
underlying assumption, the likelihood function is given by a log-normal distribution with central value J 0 bs,i and 
width <T,:, 


Ej(Ji | Jobs,P P't ) 


_ 1 _ e -(logi 0 (J;) —log 10 (Jobs,,)) 2 /2<rf _ 

Inf 10)./ obSji \/‘2na, 


(5) 


We note that Ackermann et al. [10, 13] used a different form for the J-factor likelihood function with Cj(Ji | J 0 bs,i> 0») = 
'P(Ji) — i.e., a log-normal posterior. The J-factor likelihood function used in this work differs in the substitution 
of nominal J-factor (J 0 bs,i) for the true J-factor (R) in the denominator of Equation (5). The log-normal likelihood 
formulation has several advantages over the log-normal posterior used in Ackermann et al. [10, 13]. When interpreted 
as a sampling distribution for J 0 bs,i, it is properly normalized for all values of R. The maximum likelihood estimator 
Ji for the J-factor also coincides with its nominal value J 0 bs,i from the stellar kinematic analysis. 
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FIG. 4. Comparison of the 95% CL upper limit on (ov) with its true value for a set of 400 MC realizations of the combined 
dSph analysis. In each realization a DM model with jtidm = 25 GeV and a bb annihilation spectrum is injected at the locations 
of the dSphs using the set of nominal dSph J-factors. The DM cross-section is uniformly randomized in log(au) between 1CV 2 ' 
and HT 25 cm 3 s -1 . The measured J-factor (J 0 bs,i) of each dSph is randomized by sampling from the J-factor posterior. The 
color of each point indicates the square root of the DM test statistic evaluated with the same spectrum as the injected DM 
model. Points above the solid diagonal line represent realizations where the upper limit covers the true value of the injected 
signal. 


To confirm that our upper limits have the correct frequentist statistical coverage we have performed a series of 
independent Monte Carlo realizations of our analysis in which we include a DM signal. In these simulations the 
true J-factors are fixed to their nominal values while the measured J-factors are randomized by sampling from a 
log-normal approximation to the J-factor posterior of each dSph. Figure 4 shows the upper limits on (av) from one 
set of realizations simulated with a bb annihilation spectrum and TOdm = 25 GeV. Under the assumption that the 
J-factor posterior is a good representation of the sampling distribution for J-factor measurements, we find that our 
statistical methodology produces the correct statistical coverage for a 95% CL upper limit. 


Hybrid Bayesian Analysis 


The main results of this work are evaluated with the delta-log-likelihood method, a fully frequentist statistical 
approach. In constructing the likelihood function of the delta-log-likelihood analysis, a central assumption is that 
the posterior distribution function is a good approximation to the J-factor sampling distribution. If this assumption 
holds, then our limits have the correct frequentist statistical coverage. To determine the robustness of our results 
to the choice of statistical methodology, we have performed an alternative analysis based on a Bayesian statistical 
approach in which we marginalize over the posterior distributions of the J-factors. For each target we use the same 
LAT likelihood function as for the primary analysis but set Cj = V(Ji). We then marginalize over the J-factors to 
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derive a one-dimensional marginal posterior density, 

V((fjv)] = SUi Aj((w),w DM , Ji)ir((<jy))dJi 

I n,: A,:((cre) , ) m D m, J-)n((av)')d(crv)'dJ( ’ 


( 6 ) 


where ir((av)) is the prior for (crv) and 


\i{(<7v),mnM,Ji) = Ci((crv),mDM,Ji,aLi) x 'P(Jj) (7) 

is the product of the likelihood for target i with its J-factor posterior. Given the marginal posterior density of 
Equation (6), we derive an upper limit by finding the value (crv )o that satisfies f^ av ^ 0 V((av))d(av) = p where we use 
p = 0.05 to define a Bayesian equivalent to the frequentist 95% CL upper limit. 

An important consideration for the Bayesian analysis is the choice of the prior distribution, n((av)), which is needed 
to evaluate the posterior density in Equation (6). In order to choose a prior that minimally influences our inference on 
(crv), we consider the class of non-informative priors derived according to Jeffreys’ rule [54]. As two approximations 
to the Jeffreys’ prior for our likelihood, we take the Jeffreys’ prior for the mean of a Gaussian distribution of known 
width, the uniform prior with Tc(fi) = 1, and the Jeffreys prior for a Poisson distribution, 7r(/x) = /r -1 / 2 , which we refer 
to here as the Poisson prior. The uniform prior should be applicable when the expected background is large relative to 
the signal and the LAT sensitivity is background-limited. In this regime the likelihood function of the LAT data given 
the model asymptotically approaches a Gaussian distribution. On the other hand, the Poisson prior is applicable when 
the expected background is negligible and the likelihood is well approximated by a Poisson distribution. For spectral 
models of WIMP annihilation through quark or lepton channels, the background- and signal-limited sensitivity regimes 
correspond to models of low and high mass, respectively. 

Figure 5 compares limits for the bb channel calculated with the delta-log-likelihood and Bayesian analyses. In 
this comparison we calculate two sets of Bayesian upper limits using the (crv) marginal posterior (Equation 6) and 
substituting the uniform and Poisson priors for tt((ctv)). We find that the Bayesian upper limits are in good agreement 
with the limits of the delta-log-likelihood analysis when the appropriate prior is chosen for the form of the likelihood 
on (crv). For DM masses below 100 GeV where the likelihood is well approximated by a Gaussian, the limits from 
the Bayesian analysis with a uniform prior he within 10% of those from the delta-log-likelihood analysis. At higher 
DM masses, a similar level of agreement (10-20%) is observed when comparing the delta-log-likelihood limits to the 
limits evaluated with the Poisson prior. We note that these changes are comparable to or smaller than the effect of 
the systematic uncertainties considered in the following sections. We conclude that our upper limits are robust to the 
choice of statistical methodology used to model the J-factor uncertainties. 


Systematic Uncertainties 

The dominant systematic uncertainties of this analysis arise from incomplete knowledge in three areas: the LAT 
instrument response, the Galactic diffuse gamnra-ray background, and the distribution of DM in the dSphs. To 
estimate the impact of these uncertainties, we repeat our DM search using varying assumptions intended to encompass 
the range of possibilities in each of these areas. Below we address systematics associated with the IRFs and diffuse 
background model, which both affect constraints at the 10% level, with the latter becoming less relevant for hard DM 
spectra (todm > 100 GeV). Systematics associated with the J-factors are addressed in the following section, while 
here we quote the maximum deviation from our fiducial NFW model, which occurs when assuming a cored Burkert 
density profile [55], 


Pdm(0 


Pdr s 


(r s + r)(r 2 + r 2 ) ‘ 


( 8 ) 


The J-factor systematic uncertainty has a greater impact than that of the IRF or diffuse models, approximately 35% 
at 100 GeV. We provide a summary of the systematic uncertainty as a function of DM mass and annihilation channel 
in Table II. 

In addition to the standard model of interstellar gamma-ray emission for the LAT, we consider eight alternative 
models to sample a fairly wide range of possibilities for the diffuse gamma-ray background [56]. Although we can 
vary parameters within our background models, there are no doubt sources of gamma-ray emission that remain 
unmodeled. It was observed by Ackermann et al. [13] that the TS distribution from random blank sky locations 
deviated from statistical expectations, suggesting an incomplete background model. This indicated that a rescaling 
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FIG. 5. Comparison of upper limits (bb channel) for the combined analysis of 15 dSphs as derived with the delta-log-likelihood 
analysis (solid line) and the Bayesian analysis performed with a uniform (dashed line) and Poisson (dot-dashed line) prior. The 
lower panel shows the ratio of these curves to the limits for the delta-log-likelihood analysis. 


was necessary when converting from TS to significance, effectively lowering the sensitivity of the study. One large 
class of objects known to be unmodeled are sub-threshold point sources, i.e., those which contribute gamma rays but 
are not significant enough individually to be included in a catalog. It has been speculated that these give rise to the 
larger than expected rate of type I errors (false positives) that skew the TS distribution relative to the expectation 
from Poisson statistics [13, 57]. 

Figure 6 shows the distribution of TS obtained from the analysis of randomized ROIs when the data are analyzed 
using the 2FGL and 3FGL catalogs. We additionally analyze simulated ROIs with the 3FGL catalog using an 
input model for the simulations that includes 3FGL sources and our templates for the Galactic and isotropic diffuse 
backgrounds. Using the 3FGL, which roughly doubles the number of modeled sources, brings the TS distribution 
closer to the asymptotic expectation. However, a significant deviation with respect to the asymptotic expectation 
from Chernoff’s theorem is still observed, indicating that additional unmodeled components may still be present in 
the data. 

The uncertainty in the LAT response is bracketed by using IRFs that are maximally and minimally sensitive to 
our signal. The maximally sensitive set has a greater effective area, narrower PSF, and accounts for dispersion in the 
reconstructed energy. The minimally sensitive IRFs are the opposite. The effective area is set at the boundaries of 
the envelope described in Ackermann et al. [49], while the energy dispersion and PSF width are scaled by ±5% and 
±15%, respectively. 


J-factor Uncertainties 

We have described a statistical methodology to account for uncertainties on the J-factors when deriving limits 
on the DM annihilation cross section. This procedure captures statistical uncertainties on the J-factor arising from 
the analysis of stellar velocity dispersion data. As implemented in the likelihood, it is not intended to account for 
additional systematic uncertainties. Such systematic uncertainties on the J-factor include parameterization of the 
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TABLE II. Effect of systematic uncertainties for various WIMP masses and channels reported as a symmetrical relative deviation 
from the combined 95% CL upper limits. 




10 GeV 

100 GeV 

1 TeV 

10 TeV 


IRFs 

6% 

10% 

n% 

n% 

e"*"e 

Diffuse 

12% 

6% 

3% 

2% 


J-factor 

29% 

31% 

17% 

16% 


IRFs 

6% 

10% 

11% 

11% 

+ 
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Diffuse 

13% 

6% 

3% 

2% 


J-factor 

28% 
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IRFs 

7% 

9% 
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6% 
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J-factor 
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6% 
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TS TS 


FIG. 6. Cumulative distribution of TS from 7500 randomized blank-sky regions fit with a DM particle annihilating through the 
bb channel. Shaded bands indicate the one sigma uncertainties on the cumulative fraction, which are highly correlated between 
bins. Left panel: Distributions for a DM mass of 25 GeV evaluated from Monte Carlo simulations (black line) and LAT data 
analyzed with the 3FGL (blue line) and 2FGL (red line) point-source catalogs. Right panel: Distributions for DM masses of 
25 GeV, 100 GeV, and 1000 GeV evaluated with LAT data analyzed with the 3FGL point-source catalog. As illustrated by this 
plot, the TS distribution depends on the parameters of the WIMP spectral model (annihilation channel and mass). Harder 
spectral models (e.g., those with higher mass) have a distribution that lies closer to the asymptotic expectation from Chernoff’s 
theorem. 


DM profile and the choice of priors for the profile parameters. While previous studies have shown that the derived 
J-factors are robust against these systematic uncertainties for dSphs with large stellar data sets [58], it is nonetheless 
important to quantify their impact. To assess the impact of systematic uncertainties in the J-factor derivations, we 
examine a set of four alternative J-factors derived by various fitting methods. 

The first set of alternative J-factors comes from the recent analysis of Geringer-Sameth et al. [9] assuming a 
generalized NFW profile with non-informative priors on its parameters. We also examine the J-factors derived by 
Charbonnier et al. [61] assuming a generalized Hernquist profile with uniform priors. Additionally, we perform our own 
alternative analysis following the procedure of Essig et al. [59] assuming a simple NFW profile with non-informative 
priors on the scale radius and scale density, and a velocity anisotropy parameter that is assumed to be constant with 
radius. Lastly, we show results derived from the multi-level modeling approach of Martinez [8] assuming a cored 
Burkert profile as presented in Ackermann et al. [13]. 
















16 


For each of these alternative sets of J-factors, we re-derive the limit on the DM annihilation cross section in the 
context of the LAT data. For cuspy spatial profiles, the spatial template of the DM distribution has little impact on the 
LAT analysis. Thus, for the first three sets of alternative J-factors we only alter the nominal J-factor and associated 
uncertainty. When assuming a Burkert profile, we use the full spatial profile of the assumed DM distribution (the 
change in spatial profile affects the limits by < 5%). Since the analyses of Charbonnier et al. [61] and Geringer-Sametli 
et al. [9] do not include all of the dwarfs used in our analysis, when a dSph is missing from one of these data sets we 
assign it the nominal J-factor and uncertainty from Table I in the main text. When asymmetric errors are given for 
the best-fit J-factor, we use the geometric mean to set the width of the log-normal J-factor likelihood. 

The resulting change in the upper limit on the DM cross section is shown in Figure 7. The mass dependence of the 
curves in Figure 7 reflects the fact that by changing the J-factors we change the relative importance of each dSph, 
leading to an interplay between the LAT data and the assumed J-factors. Unsurprisingly, the largest change in the 
upper limit comes from requiring a cored Burkert profile. This increases the upper limit by a factor of 20-40% with 
respect to the nominal limit (this is slightly larger than was observed by Ackermann et al. [13]) and is what we quote 
in Table II as the overall J-factor systematic uncertainty. The J-factors derived by Charbonnier et al. [61] and the 
alternative analysis with non-informative priors both yield slightly smaller changes in the limit. Finally, we observe 
that the J-factors from Geringer-Sameth et al. [9] are most similar to the nominal J-factors and result in differences 
of 5-10%. 

The combined limits presented here include both classical and ultra-faint dSphs. Bayesian hierarchical modeling 
sets rather tight constraints on the J-factors of the ultra-faint dSphs as members of the dSph population; however, 
stellar kinematic data yield larger uncertainties on ultra-faint dSphs when analyzed individually. To assess the 
maximum impact of mis-modeling the ultra-faint dSphs, we split the dSph population into ultra-faint (Bootes I, 
Canes Venatici II, Coma Berenices, Hercules, Leo IV, Segue 1, Ursa Major II, Willman 1) and classical (Carina, 
Draco, Fornax, Leo II, Sculptor, Sextans, Ursa Minor) galaxies. For soft annihilation spectra (e.g., the bb channel for 
DM with mass < 100 GeV), the classical and ultra-faint populations yield comparable limits, each ~40% worse than 
the combined limit. For harder annihilation spectra with spectral energy distributions that peak above 10 GeV, the 
limits from the ultra-faint population are roughly comparable to the combined limits, while the classical dSphs yield 
limits up to five times weaker. Considering only the classical dSphs, models with the thermal relic cross section are 
excluded for slightly lower masses (<80GeV). 


Annihilation Channels 

WIMPs may annihilate through a variety of Standard Model channels. For the quark and boson channels, the 
resulting gamma-ray spectra are all similar and largely depend on todm- The three leptonic channels have harder 
spectral energy distributions with a peak in energy flux that is closer to todm- We perform our analysis for six 
representative annihilation channels (bb, r+r - , e + e~, W + W~, and uu) and for each we assume a 100% 

branching fraction. The resulting constraints, shown in Figure 8, are similar to the bb and t + t~ channels depicted in 
the main body of this work, except for the e + e _ and /i + /r” channels which are somewhat higher. 
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FIG. 7. Change in the limits derived for the DM annihilation cross section under the assumption of alternative sets of J-factors. 
Alternative J-factors are taken from Geringer-Sameth et al. [9] and Charbonnier ef al. [61]. Non-informative priors are used 
to derive J-factors following the procedure of Essig et al. [59]. Burkert J-factors are derived using the multi-level modeling 
approach of Martinez [8] and are taken from Ackermann et al. [13]. 
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FIG. 8. DM annihilation cross-section constraints derived from the combined 15-dSph analysis for various channels. 














